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A typical problem with Monte Carlo simulations in statistical physics is that they do not allow for 
a direct calculation of the free energy. For systems at criticality, this means that one cannot calculate 
the central charge in a Monte Carlo simulation. We present a novel finite size scaling technique for 
two-dimensional systems on a geometry of LxM, and focus on the scaling behavior in M/L. We show 
that the finite size scaling behavior of the stress tensor, the operator that governs the anisotropy of 
the system, allows for a determination of the central charge and critical exponents. The expectation 
value of the stress tensor can be calculated using Monte Carlo simulations. Unexpectedly, it turns 
out that the stress tensor is remarkably insensitive for critical slowing down, rendering it an easy 
quantity to simulate. We test the method for the Ising model (with central charge c = |), the 
Ashkin- Teller model (c = 1), and the F-model (also c = 1). 
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1 ■ I. INTRODUCTION 

There exist basically two methods to obtain numerical information on two-dimensional critical systems. In the 
, transfer matrix method one calculates the largest eigenvalue of the transfer matrix and thus one finds the free energy 
of the system on a L x oo cylinder. From the theory of conformal invariance one knows that this free energy is 
Q | related to the central charge c as 

O , By introducing appropriate seams on the cylinder, that alter the cyclic boundary conditions, one also obtains some of 
the leading critical dimensions x. Here one uses the result from conformal invariance that the central charge c from 
T — I ' a system with a seam is given by 

> '■ 

00 ' c= c - 12a: . (2) 

OS : 

i After the first paper 3] that used this technique this method has become very popular. 

The advantage of the method is its high numerical accuracy. A distinct disadvantage is that the method is limited 
to rather small values of L, since the required storage capacity and computer time increase exponentially with L. 
This also implies that the method is limited to discrete spin systems. Both limitations are lifted but exchanged for a 
— «, | loss in numerical accuracy in the Monte Carlo transfer matrix method A promising |||(| method to study the 

■ transfer matrix for discrete spin systems for large L is density matrix renormalization j?j . 

The second method that has extensively been used to obtain numerical information on two-dimensional critical 
I . systems is the standard |1 Monte Carlo (MC) method. It can be used for fairly large (L x L) system sizes, both for 
discrete and continuous spin systems. Critical exponents can be extracted from the finite size scaling behavior of the 
fluctuations in critical quantities like energy and order parameter. 

However, since the free energy cannot directly be measured in a MC simulation, there exists at the moment no MC 
method to evaluate the central charge. This is unfortunate since this quantity plays such an important role in the 
determination of the universality class of two-dimensional systems. 

In this paper we want to complete this palette of existing numerical methods by presenting a direct MC method to 
evaluate the central charge. We do this by constructing an operator on the lattice that in the scaling limit represents 
the stress tensor T. The stress tensor is an operator that is connected with the anisotropy of the system; when one 
allows for anisotropy in critical models, critical points in the phase diagram become critical lines. The whole of such a 
critical line falls into the same universality class, and movements along the line are governed by a marginal operator, 
having its critical dimension x = 2. This anisotropy operator is the stress tensor T, and can be defined for any critical 
model. It is, in the language of conformal invariance, the second descendant of the identity operator. The expectation 
value of T on a L x M torus is known from conformal theory and contains in particular the central charge c. By 
comparing our MC results, as a function of M/L, with this formula we obtain the central charge and the leading 
critical dimensions. 
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An additional advantage of the method turns out to be that the autocorrelation function of T in the MC sequence 
decays much faster than that of a critical quantity like the energy. This can only partly be attributed to the fact that 
T is not a relevant but only a marginal operator. This means that there is no need to invoke more sophisticated MC 
methods like that of Swendsen-Wang M to avoid critical slowing down. 



II. CONFORMAL INVARIANCE OF CRITICAL FIELD THEORIES 



Besides being invariant against a rescaling of the length parameters, critical models are believed to be conformaly 
invariant as well: their large scale behavior is invariant against transformations that correspond locally to a rotation 
and a rescaling. Such transformations are called conformal transformations. From this symmetry, present at criticality, 
follows the structure of the Hilbert space for a large part, at least in the case of two-dimensional models. We will 
summarize some results that we need in the sequel; more details can be found in the review by Cardy fil. 




FIG. 1. The torus geometry on which the conformal field theory is defined. Dimensions of the torus are L x M, the boundary 
conditions are such that the indicated points are identified: they are cyclic in the horizontal direction and cyclic with a shift 
over M x in the vertical direction. 

In this section, we will be concerned with a system defined on a 'skew' torus; its dimensions are L x M, and 
boundary conditions are cyclic in the horizontal direction and cyclic after a shift over M x in the vertical direction, as 
in Fig. [l| Denoting the transfer matrix of the system with exp(— H), where H is the Hamilton operator, its partition 
function on such a geometry is 

Z = ^(j\e- MH e^ p \j), (3) 
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where the summation is over all configurations |j) in a row. The states \j) make up the Hilbert space on which 
the transfer matrix (or operator) exp(— H) acts. P is the momentum operator, the generator of translations in the 
horizontal direction. The Hamilton operator H of the model is the generator of translations in the vertical direction 
and commutes with P. When the transfer matrix is suitably defined, there exists an orthonormal basis of the Hilbert 
space, consisting of eigenstates of the Hamiltonian and of the momentum operator. From the theory of conformal 
invariance follows that these eigenstates with their eigenvalues are closely related to the critical dimensions of the 
model. 

There turns out to be a set of fundamental operators present in the theory, that are indicated as L n and L n for 
n e Z. They satisfy the celebrated Virasoro algebra. The Hamiltonian H and the momentum operator P can be 
expressed in terms of Lq and Lq as follows, 

27T — 7TC 

H = E L+—(L + L )-—, 

L 6L (4) 

P=fl(L -L a ). 
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Here c is the central charge of the model, and Eq is the bulk groundstate energy of the Hamiltonian, which is usually 
subtracted from it. The eigenstates of the Hamiltonian and of the momentum operator are labeled as |A + m, A + m), 
with the relations 



Hence 



L |A + m, A + fh) = (A + m)|A + m, A + fh) 
Z |A + m, A + 771) = (A + m)\A + m, A + fh) 

2vr / . c 



(5) 
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(iJ-S )|A + m,A + m) = — (A + A + m + fh - — J | A + m, A + m), (6) 

- 2tt - 

P\A + m,A + fh) = — (A - A + m - fh) |A + m,A + m). (7) 

L 

The states | A, A) with m = fh = are called primary states and the states with m and/or fh unequal to zero are 
their conformal followers. The values of A and A are related to the critical dimensions x and spin indices I of the 
operators of the theory according to 

x = A + A + m + m, 

(8) 

I — A — A + m — to. 

The appearing values of A and A from the primary states, together with their multiplicity (the level of their de- 
generacy) as well as the multiplicities of their conformal followers, determine the full structure of the Hilbert space. 
The values of the critical dimensions and their multiplicities are universal. This implies that the partition function 
(divided by its bulk value that results from Eq) considered as a function of M/L, is universal in the scaling limit of 
L and M lar ge. 

From Eq. ( |ll| ) and Eq. (||) follows that the diagonal element of exp(— MH) exp(iM x P) for the state |A + m, A + fh) 

is 

Q-c/24g-c/24QA+mgA+m (g) 

with 

/ 2-kM 2mM x \ , . 

Q = exp f — + — ^ J (10) 

and Q the complex conjugate of Q. Summing over all diagonal elements yields the partition function of the model. 
Let us label the critical dimensions with j, then 

Z/Z hulk = Q- c ' 2A Q- c / 2i NjQ Aj+mj Q Aj+ ™ 3 . (11) 

3 

This expression is called the universal expression for the partition function, and contains the central charge c, the 
values of the critical dimensions Aj + rnj and Aj + fhj as well as their multiplicities Nj. In the limit M/L — > oo, it 
yields the well known finite size scaling relation for the central charge, used in transfer matrix calculations, 

f(L)=f(oo)-^ (12) 

where 

M/L— >oo ML 

The free energy, however, is not directly accessible in MC simulations. Consequently there exist no MC results that 
yield the central charge of critical models. The stress tensor, however, is an operator that is closely related to the 
Hamilton operator and it this operator (or rather its expression in terms of spin variables) that actually is accessible 
in MC simulations, in contrast to the free energy. We will show this in the sequel. 

Conformal invariance in critical field theories states that the action (or in statistical mechanics terms: the classical 
interaction) is invariant against conformal transformations. The change in the action for non-conformal transforma- 
tions is determined by the stress tensor T(r), 

T(r) — (Txx{r) T xy (r)\ , 



3 



where T(r) is a symmetric, traceless tensor. Hence T xx (r) = —T yy (r) and T xy (r) = T yx (r). Usually, one defines the 
independent components of T as 

T ( u > v ) = ^ [Txx(u, v) - iT xy (u, v)] , 

I (15) 
T{u, v ) = 2 w) + iT xy (u, v)] ■ 

Here u and w are the position coordinates on the torus. The dimension (A, A) of T and T are (2,0) and (0,2) 
respectively. Hence their critical dimension x — 2 and their spin indices are I = ±2. So the stress tensor is a marginal 
operator. Its components T and T can be expressed in terms of the fundamental Virasoro operators L n and L n as 
follows: 

2n\ c 
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This expression is valid on the torus geometry; u (v) is the horizontal (vertical) position on the torus. The Virasoro 
operators L n and L n with n^O play the role of raising and lowering operators for the states |A + m, A+m). Because 
these states are orthonormal, only L and L have nonvanishing contributions in the expression for the expectation 
value of T. From the expression ( p7j| ) and the eigenvalue equations (|^), the expression for the expectation value of the 
stress tensor can easily be calculated. The expectation value is 



(T) = \ Y, + m ^ h i + %l T \^ + m ^ A o + %>• (17) 
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Substituting Eq. (|l6| ) and the expression for the partition function (^JJ) yields 



7 = I — ) * (- ^' Nj {A] + m ° ] Q A '" Q ' 



Ly \ 24 y A '" (; A m / 

(18) 

. , , . i ™ A n-A,-m<n-A,-mi \ 



2n\ 2 / c Z) ■ -Nj (A J a ™ A n A "' ° A "' 



(T) - , 

i ' ' '24 V XjQ A - '" Q A '" 

We will be mostly concerned with the diagonal elements T xx = ~T yy of the stress tensor. The expression for T xx is 

■2- Y ( c Ej Nj (Aj + mj + Aj + %) Q A Q A 



(Txx) (^ L J ^ 12 V .Y,<7 A »• Q A J ^ 

This expectation value can be written as the derivative of the free energy with respect to the aspect ratio M/L, as 

Upon taking the derivative, the volume ML of the system is kept constant. Like magnetization and magnetic field, 
the 'field' M/L is the external field conjugate to the operator that is the stress tensor. Therefore, (T xx ) couples to 
the anisotropy of the system. 



III. A LATTICE REPRESENTATION OF THE STRESS TENSOR 



In conformal field theory, the stress tensor is an operator that is quite abstract. It is defined only after the lattice 
model has gotten its continuum limit. For lattice models, however, the stress tensor can easily be defined as well. This 
lattice representation of the stress tensor must thus have the scaling behavior predicted by expression ( p~9[ ) . Below we 
will illustrate the construction of the lattice representation of the stress tensor for the Ising model, but first we give 
a more general way to proceed. 
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A. Constructing the stress tensor 



Construct, for a lattice model, an operator t(r) as an expression in the local, fluctuating field(s), such that: (i) 
t(r) transforms as a second rank tensor (in particular, i(r) picks up a minus sign under a rotation over 90°); and (ii) 
t(r) has the same symmetry as the interaction energy of the model under study. In general, this means that t(r) is 
invariant under global spin flips or spin rotations. 

If one now expresses i(r) in terms of scaling operators it is clear that the operators that occur should all be tensors 
that change sign under rotations over n/2, i.e., they have / = A — A = ±2, ±6,.... Since A and A are always 
non-negative, it follows that the scaling dimension x of the appearing scaling operators all have x > 2. The marginal 
case, having x — 2 and I — ±2, is in fact the stress tensor, all other operators in the expansion are irrelevant. To be 
more precise: in this general case, the operator t(r) couples to both independent components T xx (r) and T xy (r) of 
the stress tensor. As will become clear below, however, t(r) can easily be defined such that it couples to T xx (r) only. 
In that case, one has 

f(r) = a T xx (r) + ■ ■ ■ , (21) 

where the dots represent irrelevant operators. The requirement (ii) guarantees that i(r) and T xx (r) share the same 
interaction symmetry, so that the coefficient a does not vanish by symmetry. Constructing operators t(r) can, as we 



shall see, be done in several ways, but all choices yield an expansion (21), albeit with different values of a. 

Having constructed the operator t(r) one can evaluate its average in a MC simulation on a geometry of L x M, for 
several values of M/L and L large. The result should follow the universal expression for T xx (r) as 

(t(r))=a{T xx (r))+0(±j), (22) 

where the expression for (T xx (r)) given in Eq. ( |l^) is proportional to 1/L 2 and dominates the second term that has 
uj > 2. Hence we can fit the M/L dependence of the left hand side against Eq. ([l9]), obtaining, in particular, the 
central charge c. 



B. The stress tensor for the Ising model 

We will illustrate the construction of the discrete stress tensor t(r) in the case of the Ising model. The starting 
point is the close connection between stress tensor and anisotropy. Let us therefore start with the anisotropic action 
A of the ordinary, square lattice Ising model, 



A — — ^ (j x SijSi+xj + J y Si t jSi t j+i*j , (23) 



ij 



where the couplings (J x ,J y ) allow for anisotropy The isotropic critical point is J x = J y = J c = ^ln(l + v2), but 
this point becomes a critical line when unequal values of J x and J y are allowed for. 

The central notion here is that in the scaling limit, anisotropy amounts to a rescaling of the length parameters x 
and y with a different scaling factor. Hence, in the scaling limit, the anisotropic model with (J x , J y ) behaves as the 
isotropic model with rescaled length parameters x and y, 



y' = e+V 



(24) 



The value of A in this equation determines the values of J x and J y . In this way, Eq. (|24|) fixes the parameteri- 
zation [J X (X), J y (X)] of the critical line with A. The isotropic point has A = with J x (0) — J y (0) — J Cl and the 
parameterization obeys J X (X) = Jy(— A). 

On a finite geometry L x M, this anisotropic rescaling means that the volume ML of the system remains untouched, 
but that the aspect parameter M/L scales according to 
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In the scaling limit, therefore, the partition function with the anisotropic action of equation (|23|), which we call 
Z(X,M/L) depending on A, equals that of the isotropic Hamiltonian with a rescaled aspect ratio M' / L' , 

ZU^).Z(,.0.^). (2(5) 



L J V L . 

A general movement in the phase diagram is performed by a scaling operator. A renormalization transformation 
is isotropic, which implies that there can be no renormalization flow along the critical line (J x ,J y ). This implies 
that the scaling operator that governs the movement along this line must be invariant against a renormalization 
transformation, i.e., it is a marginal operator having its critical dimension x = 2. 

In the case of the Ising model, the action (|2^) immediately shows which operator this must be. Write the action as 
a symmetric part plus a part that determines the anisotropy, 

A = A c — ^~^(Jg(A) — J c ) SijSi+ij + (J y (X) — J c ) SijSij+i, (27) 

ij 

where A c is the action at the isotropic critical point J x = J y = J c . Expanding up to first order in A, this expression 
can be written as 

A = A C + X^t xx {i,j), (28) 

ij 

where t xx (i,j) is the lattice representation of the stress tensor, 

txx(i,j) = -JL{ty(SidSi+i,3 ~ S%,jSi,j+ij. (29) 

Here we used the symmetry property J X (X) = J v (—X). The operator t xx (i,j) governs the anisotropy of the system. 
This lattice representation of the stress tensor ( p9[ ) for the Ising model was already known for a long time JlO| ; the 
value of J' x {0) = \y/2. 

In fact, the operator in Eq. ( p9| ) is written as t xx because it is one of the two components of the full stress tensor 
t a p(i,j). This t a p(i,j) has the same properties as the field theoretical stress tensor: it is a second-rank, symmetric 
traceless tensor. The other component, t xy (i,j) can be written as 

t X y(i,j) — ~J X y(^SijSi+l,j + l — Si t jSi+lj — lj, (30) 

with a certain prefactor J xy , which will be different from J x (0), because it couples next-nearest neighbor spins instead 
of nearest neighbors. The off-diagonal elements of the discrete stress tensor couple to the anisotropy in the diagonal 
directions. 

It is this operator t a/ 3(r) that appears in the previous subsection. It is constructed such that it behaves as a 
second-rank symmetric tensor with the same symmetry under global spin flips as the interaction energy itself. Of 
course, this version of t a p(i,j) is not the only possible one; it can also be defined with further neighbor interactions. 

The precise connection between the discrete variant t(r) of the stress tensor and its field theoretical counterpart is 



obtained by taking the derivative of Eq. (E6[) with respect to A at A = 0. Using Eq. (20), this yields 



J' x {0)J2(SijS l+hj - S itj S ilj+l } - — (T xx (u,v)), (31) 



where (T xx (u, v)) is the expression ( |l9| ) for the expectation value of the diagonal component of the stress tensor. Note 
that this expression is a universal function of central charge, critical dimensions, and their multiplicities. The value 
of J x (0), however, is in general unknown, such that we will have to include it as a fit parameter. 

Expression ( |3l| ) combined with Eq. (|l9|) yields the relation that is central to this work: it expresses the expectation 
value of the lattice representation of the stress tensor in the universal quantities that we want to know. As we will 
use a rectangular geometry, without the shift in boundary conditions, we put M x = and obtain 



(SijSi+ij - SijS irj+ i) — a ( — ) — ^ 3 — 7 - M L I , (32) 
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FIG. 2. The expectation values of stress tensor (1) for the Ising model, defined in Eq. (ffl|), as a function of the aspect 
parameter M/L. From high to low, the plots show the expectation values from system dimensions L running from 4 to 10. 
The lines are the result of the fit against Eq. (^) together with a correction to scaling term (^jj). Note that, for each value of 
L, the stress tensor at M/L = 1 is zero by symmetry. 

where 1/cv = it J^(0). See Fig. |^ for an example of the functional dependence. The prefactor a is the same a 
that appears in Eq. (21). The physical interpretation of it is given by the relation 1/a = irJ x (0); it determines the 
'amount of anisotropy' that the system obtains, once the stress tensor is switched on. Note that a is non-universal; it 
depends on the precise definition of the model, as well as on the definition of the stress tensor. The other quantities 
present in Eq. (|32|), however, are the central charge, the critical dimensions, and their multiplicities, and those are 
all universal. There is an infinite number of critical dimensions, but only a limited number of these is 'small', say, 
less than two. For large enough values of the aspect ratio M/L only a limited number of critical dimensions have 
nonvanishing contributions to Eq. (§3) such that a fit of this expression against MC data must be feasible. Typically, 
we will take M/L > 1. 



IV. FITTING THE MONTE CARLO RESULTS 

Fitting the MC results against the universal expression for the stress tensor (^2|) requires a decent fit program, as 
the number of fit parameters is quite large, and requires some theoretical reflection on the model as well. We will 
deal with the use of the universal expression (f32|) and the corrections to scaling in different subsections. 



A. The universal expression for the stress tensor 

The number of critical dimensions Xj that appears in the universal expression (^2|) for the stress tensor is infinite, 
which clearly is an infeasible number of fit parameters. Most of the dimensions, however, are large. Their contribution 
to Eq. ( pS2] ) goes as exp(— 2ttxM/L), so if we limit the calculations to values of the aspect ratio M/L that are not 
too small, most of the dimensions Xj have a vanishing contribution. Performing some preliminary MC simulations 
suggests a reasonable lower bound to M/L. Typically, we took M/L > 1. The upper bound on M/L is determined 
by the value of M/L where (t xx (r)} reaches its asymptotic value. To determine a reasonable upper bound on M/L, 
the same preliminary simulations can be used. The asymptotic value of (t xx (r)) is 

(Mr)>| M/i _ = a(f) Y2 (33) 

This expression also shows why the asymptotic value itself is not sufficient for the determination of the central charge: 
it only gives an estimate of ac instead of c. 

Performing simulations to obtain the expectation value of the stress tensor between these bounds on M/L is in 
principle sufficient to extract the desired quantities by fitting the results against expression (|32]). Typically, we take 
three or four critical dimensions into account, the identity dimension xq = 0, present in any critical model, and two 
or three nontrivial ones Xj. For these dimensions, the multiplicities Nj must be specified in the expression as well. 
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A multiplicity is always integer, and when it is larger than one, the corresponding critical dimension x is degenerate, 
and thus an additional symmetry is present in the model. Some theoretical reflection on the model often is sufficient 
to reveal such symmetries. Another possibility is to perform fits with different values of the multiplicities and to 
choose the set that gives the best fit. 

The lowest appearing values of the dimensions x — A + A correspond to primary fields. As noted in Sec. to each 
primary field belongs a tower of conformal followers or descendants, that have values of the critical dimensions that 
differ by an integer from that of the primary field; they are A + m and A + fh with m, m G N. The first descendant 
of a scalar primary field 0(r) with dimensions (A, A) is VO(r), a vector field that has two components; one having 
(A + 1, A) and the other having (A, A + 1). To improve the fit, we will include this first descendant into Eq. (|3^). 
This inclusion introduces no new fit parameters; the value of its critical dimension is x + 1 when x is the critical 
dimension of the primary field, and this value appears twice. Hence the multiplicity of the first descendant is twice 
that of the corresponding primary field. 

Having fixed the multiplicities in expression (|32[), this leaves us with four or five free parameters: the prefactor a, 
the central charge c, and two or three nontrivial dimensions Xj. 



B. Corrections to scaling 

The above analysis of the universal behavior of the stress tensor is valid in the scaling limit. The discrete version 



of the stress tensor t xx (r), however, is not a scaling field; as argued in Sec. Ill it can in fact be written as the 
expansion ( |2l| ) in scaling fields, of which only the first term is the true stress tensor with its universal behavior. The 
fit to this expression is treated in the previous subsection, but for smaller system sizes other terms in the expansion 
become important. The scaling behavior of these terms in L goes as with u> > 2. 

To obtain accurate results, we should include at least one of these correction terms in the expression that we fit 
against our MC results. That means that we have to perform calculations for different values of the width L of the 
system in order to be able to extract the L~ 2 behavior of the true stress tensor. 

In principle, we could proceed by performing simulations for a fixed value of M/L and increasing values of L and 
extract, by extrapolation, the part of (t xx (r)) that scales as L~ 2 . This means that, for any value of M/L, we have to 
fit 

(Ur)) = ^ + p (34) 



and have to use the values of a for each value of M/L to fit against expression (|32|). We can, however, do better. 
To this end, write the first scaling operator on the dots in expansion ( pl| ) as O(r), 

(t xx (r)) = a (T xx (r)) + (O(r)). (35) 

Now consider the general expression for the expectation value of an operator O(r) on a system with geometry L x M, 

(0(r)> = ^^<j| e-^O(r) (36) 
j 

with Z the partition function and H the Hamilton operator of Eq. (^). Using the basis | A + m, A + fh) of the Hilbert 
space yields 



where we used Eq. . Here the parameters aj(L) are the diagonal elements of the operator O(r) in the basis 
| A + m, A + fh) of the Hilbert space, 

a,j(L) = (Aj + rrij , A j + fhj\ O(r) \Aj + m^Aj + fhj). (38) 

As the basis functions | A + m, A + fh) only depend on L and not on M/L, the full M/L dependence of the expectation 
value (O(r)) is accounted for by the exponentials in Eq. (|37|). The amplitudes a,j(L) depend on L only. Taking only 
the leading correction into account, they can be written as 

"i(L) = ^, (39) 
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with the same value of ui for each of the diagonal elements [jjj. In our fit, we will only include the few most important 
critical dimensions Xj\ the identity dimensions x$ — and the first two nontrivial ones. That means that including 
the expression for (O(r)) as a correction to scaling gives only four additional fit parameters: the correction exponent 
w and three amplitudes do, a\, and ai. 

In the other case, by naively extrapolating the behavior of (t xx (r)) for large L, we have to include the two fit 
parameters u> and b for each value of M/L. The above approach thus drastically reduces the number of fit parameters. 

Still, in the complete analysis of the MC results, the number of fit parameters is quite large. Typically, we need 
four parameters from the expression (|32|), which are the prefactor a, the central charge c, and the two most relevant 



dimensions x\ and X2- For the corrections to scaling, we use the expression following from Eq. (|37|) and (39) 



«*'» = tp (40) 

giving four additional parameters, which are uj, ao, a\ and ai- 

In this way, we perform a combined fit of the MC results for all values of L and M/L, in one single fit using eight 
fit parameters. This is a large number, but the functional dependence of the formula, to be fitted for two variables 
simultaneously, is very restrictive. Especially for expression (p3), the behavior in L is restricted to L~ 2 , and the values 
of the dimensions x\ and xi appear as dimension as well as as amplitude. 

The number of critical dimensions Xj and the values of M/L that have to be included in the fit are a matter of trial 
and error. Sometimes it turned out necessary to delete some of the lowest values of M/L from the data set. Lower 
values of M/L clearly stabilize the fit, but on the other hand, including these values requires more critical dimensions 
from expression ([52j) to describe the full data set. We varied the lower bound on M/L and the number of critical 
dimensions, until the quality of the fit became high enough. 

This procedure requires a fit program that yields, apart from the values of the fit parameters and their error bars, 
also a parameter that indicates whether the fit can be trusted or not. Our case amounts to a two-dimensional fit (in 
L and M/L) using eight or ten fit parameters. The program we used is based on routines from Numerical Recipes 
jnj. The parameter that indicates the quality of the fit is called the goodness of fit Q. The value of Q lies between 
and 1, and is based on the x 2 of the fitted data. Q gives the probability that the x 2 of a certain data set exceeds 
that of the actual data set. A very low value of Q means that it is highly unlikely that the used function gives the 
correct theoretical description of the data. In our case this means that we either included values of M/L that are too 
small, or not enough critical dimensions Xj. 



V. COMPARISON WITH EXACTLY SOLVED MODELS 



In order to test the method, we performed MC simulations on some models of which the scaling behavior on 
the torus is known exactly. We chose the Ising model (with central charge c = |), the Ashkin- Teller model (with 
c = 1) and the F- model (also with c = 1). There is a line in the phase diagram of the Ashkin- Teller model that can 
be mapped, by a duality transformation and a graphical representation JT2j, exactly on the F-model. We chose to 
simulate the corresponding points in the Ashkin- Teller model and the F-model. The results, however, differ, which is 
an illustration of the importance of boundary conditions in such simulations. The duality transformation alters the 
boundary conditions, giving rise to a different behavior of both models on a finite geometry 

In case of the Ising and Ashkin- Teller models, we performed MC simulations using the standard Metropolis algo- 
rithm. For the F-model, we had to use a cluster algorithm as well (to be described below) . We performed simulations 
on a system with geometry L x M with varying values of L as well as of M/L. We sampled different versions of the 
stress tensor t xx {v), in order to obtain independent estimates of central charge and critical dimensions. 



A. The Ising model 

We carried out our simulations on the ordinary square lattice Ising model, with the action given in Eq. ( p3[ ) on 
its isotropic critical point given by J x = J y = J c = ^ ln(l + V2). The construction of the stress tensor is described 



in Sec. [II. Actually, taking different versions of the discrete stress tensor t xx (v) gives an independent check on the 
accuracy of the results. All different versions should couple to the true stress tensor T xx (r), albeit with different 
prefactors a. We chose two versions of the discrete stress tensor, one defined with nearest-neighbor couplings and the 
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other with next-next-nearest neighbor couplings: 

(1) (SijSi+ij ~ SijSij+i), 

(2) (Si,jSi+2j — SijSij+2}. 

Note that the stress tensor defined with next-nearest neighbor couplings corresponds to the off-diagonal elements of 
the stress tensor; its expectation value on the used geometry is zero by symmetry. We took the system geometry 
L X M with L varying from 4 to 10 and M/L varying from 1.5 to 10. 

TABLE I. Monte Carlo results for the Ising model. Stress tensors (1) and (2) refer to the definition in Eq. (ptl|). Values of 
the prefactor a, central charge c and the first two critical dimensions xi and X2 are given and compared with their exact values. 
Errors in the last digit are given between parentheses, iv is the power of the l/L correction, and g.o.f. is the 'goodness of fit'. 
In case of stress tensor (1), the prefactor a is known exactly JfOl. 





Stress tensor (1) 


Stress tensor (2) 


Exact 


a 


0.450 (2) 


1.277 (3) 


\/2/2 = 0.4501. .. a 


c 


0.500 (2) 


0.498 (1) 


1/2 


Xi 


0.1254 (6) 


0.1256 (4) 


1/8 


X2 


1.0 (4) 


1.1 (4) 


1 


UJ 


4.3 (1) 


4.29 (8) 




g.o.f. 


0.83 


0.97 





a only for stress tensor (1) 



The resulting expectation values were fitted against expression ( |32| ) together with a correction to scaling term of 
Eq. (f4(i|). We took two nontrivial critical dimensions x\ and X2 into account, both with multiplicity 1. The data for 
stress tensor (1), together with the results of our fit, are plotted in Fig. || to get a feeling of the behavior of the stress 
tensor. The numerical results of the fit are summarized in table ||. Even for those small system sizes, accurate results 
are obtained. This is the first determination of the central charge using a MC simulation. 



B. The Ashkin-Teller model 



A more severe test of the method is obtained by considering a model having dimensions lying closer to each other. 
The Ashkin-Teller model is a useful candidate for testing our method. It has in its phase diagram a critical line which 
can be mapped on the (exactly solved) six vertex model fl2]| . The universal partition sum of the Ashkin-Teller model 
on the torus is exactly known |I3[ ], so it can be compared with our MC results. 

The Ashkin-Teller model has two Ising spins S and P with S, P € {+1, — 1} on each lattice site, that interact with 
an action 

A = - J (SiSj + PiPj) + KSiSjPiPj, (42) 

(ij) 

where (ij) denotes a summation over nearest neighbor lattice sites. The critical line in the phase diagram that can 
be mapped on the six vertex model is parameterized by 

1 - W 

exp(2J) = 



exp(2 J + 2K) = 



W (43) 
1 + W [ ' 



l-W 



The weight W equals the Boltzmann weight of the four vertices in the six vertex model that carry a step. The other 
vertices are flat and have Boltzmann weight 1. 

The critical line of Eq. ( |43] ) is a line with central charge c = 1 and continuously varying exponents. By expressing 
the partition function of the Ashkin-Teller model in the scaling limit in terms of Coulomb gas partition functions, 
all critical exponents can be obtained. For this derivation, the reader is referred to Ref. [Q; we will only state the 
results. 

Part of the exponents varies continuously along the critical line. Their value is expressed in terms of the renormalized 
value of the Gaussian coupling g, present in the Coulomb gas partition functions. The dimensions x of the primary 
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fields are 



and g is the Gaussian coupling, given by 



2 2 

— + ^— with e, to e Z, (44) 
Zg Z 



9 = l aicsin {zw)- (45) 

The other dimensions are constant along the critical line. We chose, rather arbitrarily, the point W = 0.8 on the 
critical line for our simulations. At this point, the three most relevant dimensions are 

x\ = 0.125 (with multiplicity 2), 

x 2 = 0.2908... (with multiplicity 1), (46) 
x 3 = 0.8596... (with multiplicity 1). 

Typically, the multiplicity of the degenerate dimension x\ (which is constant along the critical line) can be guessed 
beforehand, albeit some theoretical reflection on the model is necessary. To this end, consider the expansion in scaling 
operators of S and P, 

S(r) = a s p(r) + • • • , 

P(r)= Qp? (r) + -, 1 ) 

where p(r) and q(r) are the leading (most relevant) scaling operators in the expansion. The manifest symmetry S <-> P 
of the action ([42]) implies that 

a| (KriMr 2 )> = 4 (?(ri)?(r 2 )). (48) 

Hence it follows that p(r) and q(r) share the same critical dimensions x. On the other hand, spin reversal symmetry 
S — > —S implies that 

{S{v x )P{v 2 )) = (49) 

which implies that the dominant term for |ri — r2 1 large in this expression must vanish as well. Hence 

a s a P (p(r 1 ) g (r 2 )> = 0. (50) 

This ensures that p(r) and q(r) are different scaling operators sharing the same critical dimension x. Therefore this 
magnetic critical dimension x must have multiplicity 2. Note that the second argument does not apply for energy- like 
operators like SijSi+ij, such that the energy scaling field will be non-degenerate. 

We performed MC simulations using the standard Metropolis algorithm, again on the system with geometry L x M, 
with L varying from 5 to 12 and M / L varying from 1.5 to 10. We sampled four different versions of the stress tensor 

txx (l") • 

(1) (Si.j Si+ij + PjjPj_)_xj ~ SijSij+i — PijPij+i) , 

(2) (SijSi+2,j + I ■ 1., - S h jS ij+ 2 - PijPij+2), 

(3) (SijPijSi+ijPi+ij —SijPijSij+iPij+i), 

(4) (SijPijSi+2,jPi+2,j — Si,jPi,jSi 1 j+ 2 Pi,j+2)- 

Stress tensors (1) and (2) are defined such that the symmetry between S and P spins is incorporated. 

TABLE II. Monte Carlo results for the Ashkin- Teller model, corresponding to the six vertex model with Boltzmann weight 
W = 0.8. The stress tensors (1) to (4) are defined in Eq. (J5l|) . For notation see table [l| 

Stress tensor (1) Stress tensor (2) Stress tensor (3) Stress tensor (4) Exact 

~ 0.36 (2) 0.78 (3) 0.23 (1) 0.65 (2) 

c 0.97 (6) 0.99 (3) 0.96 (4) 0.96 (3) 1 

xi 0.128 (3) 0.130 (2) 0.128 (3) 0.128 (2) 0.125 

x 2 0.33 (8) 0.34 (5) 0.34 (6) 0.35 (4) 0.2908... 

x 3 0.9 (4) 0.8 (1) 0.9 (3) 0.9 (2) 0.8596... 

cj 3.8 (2) 4.1 (1) 4.0 (2) 4.2 (1) 

g.o.f. 0.80 0.016 0.76 0.065 
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It turns out that in this case three nontrivial dimensions have to be included in the fit. This brings the total number 
of fit parameters to no less than 10. Still, relatively good results are obtained; they are summarized in table |l[ 

C. The F-model 

A nice illustration of the importance of boundary conditions is obtained when a dual version of the Ashkin- Teller 
model is considered. As stated, the critical line of the Ashkin- Teller model can be mapped exactly onto the F-model, 
using a duality transformation and a graphical representation [ O ] . On a finite system, however, this mapping affects 
the boundary conditions, such that both models with periodic boundary conditions will have a different behavior on 
the torus. 

The model we chose to consider in fact is an intermediate model between the F-model and the Ashkin- Teller model, 
and is obtained from the latter by applying a duality transformation on one of the spins S or P only. In this way, 
we obtain two coupled Ising models, defined on two interpenetrating sublattices. Both Ising models are equal; they 
interact via a nearest-neighbor coupling such that a broken Ising bond carries a Boltzmann weight W , where the 
weight W is the same as the W in Eq. ( f43| ) of the Ashkin- Teller model. The coupling between the Ising models only 
exists in the restriction that two broken Ising bonds are not allowed to cross each other. An elementary square of the 
lattice contains two spins of both sublattices; diagonally opposed spins belong to the same sublattice. The restriction 
is that at most one of the bonds over the elementary square may be broken. 

The resulting model can easily be mapped on the F-model, seen as a BCSOS height model |l4j]. To this end, the 
Ising-Bloch walls are identified with the steps, carried by the first four vertices of the F-model. To become steps, 
walls have to be equipped with an arrow; the steps have to be identified as a step up or a step down. This arrow 
assignment is simply such that two adjacent Ising-Bloch walls carry antiparallcl arrows if they belong to the same 
sublattice, and carry parallel arrows if they belong to different sublattices. 

In this way, a configuration of the two Ising models is mapped onto a configuration of the F-model, and vice versa. 
There is, however, a difference in boundary conditions on the torus. If we consider the F-model on a finite geometry 
as a height model, we have to allow for defects at the boundary. The smallest defect in the F-model is a defect of two 
unit heights, which corresponds to two steps running over the system. The corresponding Ising configuration however, 
would have one Ising-Bloch wall running over the system for each sublattice, which is not allowed when the two Ising 
models have periodic boundary conditions. Hence, for the F-model the allowed defects at the boundary are height 
differences multiples of 2, whereas in the formulation of the Ising models, the height differences at the boundary are 
multiples of 4. 

Related to these defects is a complication that arises, when one naively tries to simulate this version of the F-model 
using a single-spin Metropolis algorithm. As the updates in such an algorithm are always local, it cannot generate 
configurations with defects around the torus. The algorithm is able to generate islands of flipped spins, but such an 
island never can cross an Ising-Bloch wall of the other sublattice. This implies that the algorithm is non-ergodic; 
the part of phase space it reaches is restricted to that part that has the same defects at the boundary as the initial 
configuration. 

That does not mean that the results of the simulation make no sense. The model that results when using only the 
Metropolis algorithm is a true height model, such that on the boundaries no defects are allowed at all. This model 
renormalizes to the Gaussian model. The universal form of its partition function is known fl5|| , but behaves somewhat 
anomalously because it has a continuous spectrum of critical dimensions, that result in an integral instead of a sum 
in Eq. (|32"|). The universal partition sum of the Gaussian model is the result of this integral. Inclusion of its form in 
our fit for this model indeed yields the correct result. 

The difficulty in boundary conditions, however, can easily be overcome using a cluster algorithm, that allows for 
non-local updates of the configurations. In our simulations, we used a standard Metropolis algorithm for thermal 
equilibration, combined with a cluster algorithm ]9|,|l6|] that is able to generate defects, in order to make sure that the 
whole phase space can be reached. We performed simulations on the model with L varying from 6 to 18 and M/L 
from 2 to 5. It turned out in this case the stress tensor reaches its asymptotic value already for M/L 5. 

We sampled two possible versions of the stress tensor, 

(1) (SijS i+ 2,j - SijSij+z), . 

(2) (SijSi+sj+i — SijSi-i 

where we took into account that energy-like spin products always must couple spins of the same sublattice. The most 
simple version of the stress tensor couples nearest-neighbor spins of each sublattice, but its expectation value on the 
system geometries that we considered is zero by symmetry. Stress tensor (2) is, regarding its definition, a mix of 
txx( r ) an d ^xji(r)j but this is no problem since, on the used geometry, any t xy (r) is zero. 
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TABLE III. Monte Carlo results for the F-model with Boltzmann weight W = 0.8. Two different stress tensors are used for 
the calculation of the central charge and critical dimensions. They are defined in Eq. ( [52] ) . For notation see table |. 



Stress tensor (1) Stress tensor (2) Exact 



a 0.83 (6) 1.32 (7) 

c 1.06 (7) 1.03 (6) 1 

xi 0.291 (6) 0.289 (5) 0.2908... 

x 2 0.7 (1) 0.65 (9) 0.8596... 

u 3.2 (2) 2.8 (1) 

g.o.f. 0.28 0.27 



The fact that this model is a height model ensures that there are basically two types of operators, spin wave and 
vortex operators, with dimensions given in Eq. (|44|), that both are doubly degenerate, (cf. |l7| for further discussion.) 
Hence, the lowest critical dimensions have multiplicity 2. We fitted the resulting expectation values of the different 



stress tensors, using two non-trivial critical dimensions. The results are summarized in table III 



It is noteworthy that the prefactor a in the definition of the stress tensor is independent of the boundary conditions. 
Our fit on the simulation that only used the Metropolis algorithm (described above) yielded the same prefactors as 



those in table [II. That means that the expansion (|2l|) of the discrete stress tensor in terms of scaling fields only 



depends on local properties. 

It turned out that including values of M/L smaller than 2 destroyed the quality of the fit, yielding a far too low 
value of the 'goodness of fit'. The reason probably is that there are much more dimensions Xj present that are quite 
small and that start to become important for values of M/L smaller than 2. This can be seen from the value of X2 
that follows from the fit; it is significantly lower than the exact value of the second dimension. Apparently, in the fit 
program xi plays the role of an 'effective' dimension, incorporating the values of several dimensions in one. This casts 
doubt on the validity of the highest dimension that is given by the fit program, but is seen not to affect the values of 
the central charge c and the most relevant dimension X\ . 



VI. SIMULATION TIMES AND AUTOCORRELATIONS 



MC calculations of a marginal operator like the stress tensor typically encounter additional difficulties as compared 
to observables like energy and magnetization. The latter quantities have a relative error in MC simulations that does 
not scale with the system size, whereas this is not the case for an operator like the stress tensor; its relative error 
increases with the system size. 

This can be seen as follows: consider an operator O(r) of which we want to calculate its expectation value. Its 
scaling behavior will be dictated by a critical dimension x, 

i^<0(r))~L-*, (53) 

r 

where L is the linear system size. The error A. ( r ) m the average value is related to the number of samples N in the 
MC simulation and to the second moment of its distribution, 

^ (r) = ^7\E<°( r )°( r ')> - (0(r))<0(r')). (54) 

r,r' 

Note that N stands for the number of statistically independent MC samples. The dependence on L of the simulation 
time to reach independent samples will be discussed below. 

Typically, the double summation in the last expression has two contributions; a short range and a long range 
contribution. The short range contribution, say within a region with radius R, follows 

(O(r)O(O)) - (O(r))(O(0)) constant, (55) 

\r\<R 

for L large. The constant is roughly proportional to the radius R when it is not too large. The long range contribution, 
on the other hand, is dominated by the critical dimension x as 

(OWO(0)) - (O(r))(O(0)) ~ L 2 ~ 2x , (56) 

\r\>R 
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Now there are two cases. If the dimension x < 1 the long range contribution dominates Eq. ([54j), and the relative 
error in (0(r)) scales according to 

A °W „ _L ( 57 ) 

It is inversely proportional to the square root of the number of MC samples, but does not scale with the system size 
L. This is the usual case for, .e.g., magnetization and energy in the Ising model. In case x > 1, however, the short 
range contribution dominates the error for large L, which implies that the relative error in (O(r)) scales according to 

A °(r) ■ 1 TX-1 (za) 

We will want to obtain the same relative error for all different linear system dimensions L in our MC simulations. 
For observables having x < 1 this requires the same number of MC samples for all L. For x > 1, however, Eq. ( |5S| ) 
dictates that N ~ L 2x ~ 2 . In case of the stress tensor, having a; = 2, the number of MC samples should thus be 
proportional to L 2 . 

At first sight, it seems that this fact makes it difficult to reach large system sizes, as the simulation time is directly 
proportional to the number of required MC samples. This, however, is only partly true. The other parameter which 
determines the simulation time is the time it takes to generate statistically independent configurations. Critical 
systems are known to suffer from critical slowing down. If one uses the standard Metropolis algorithm, the typical 
time r it takes to generate statistically independent configurations increases with the system size as a power law. 

Unexpectedly, it turns out that the stress tensor is remarkably insensitive to critical slowing down. This can 
be judged from its autocorrelation function. Let us define a MC cycle as one attempted update per spin. The 
autocorrelation function of a certain observable O is defined as 

a( t ) = ( Qt ° Qt °+ f > Z ( 0t °^ f 59 ) 

The operator O is, as usual, defined as ^ r ^( r )- Here O t denotes the value of O after t time steps, where a time step 
is one cycle, i.e., one attempted update per spin. The autocorrelation function g(t) is normalized such that g(0) = 1. 
In practical situations, the number of MC cycles t between two consecutive MC samples has to be such that g(t) is 
(almost) zero. 

The observation that the stress tensor does not suffer very much from critical slowing down follows from Fig. |^. 
Here we plotted the autocorrelation functions git) for the energy and the stress tensor, in case of the Ising model at 
its critical point, for several different system dimensions. For the number of cycles t not too small, the autocorrelation 
function of the energy shows a straight line in the log-normal plot, meaning that its behavior is exponential in t. 
Indeed, the behavior of the autocorrelation functions for nearly critical systems is given by 

git) ~ exp(— t/r) for t large, (60) 

where r is the autocorrelation time. The dynamic scaling hypothesis states that the time scale r of a dynamical 
system is connected with the length scale, which is the correlation length £, and that this connection is described by 
a universal dynamic exponent z, 

r~r- (6i) 

The exponent z is believed to be connected to the dynamics of the system (in our case, by the Metropolis algorithm) 
and to be the same for all observables. For finite systems at their critical point, the correlation length £ is bounded 
by the system dimension L, such that 

t(L) ~ L z . (62) 

We extracted the values of t(L), following from the autocorrelation function of the energy in Fig. [| by fitting the 
autocorrelation functions to Eq. (|6^). For this, we removed the first data points, up to the point where the plot begins 
to show a straight line. The values of t(L) were fitted to Eq. (|62|), yielding a value for z of roughly 2. The quoted 
value in the literature faM is z sa 2.17, which is consistent with our findings. 
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FIG. 3. Plots of the autocorrelation function g(t) of Eq. (f39|), where t is the number of Monte Carlo cycles. Calculations 
are performed using the Metropolis algorithm for the Ising model at its critical point. System dimensions are indicated in the 
figure, (a) Autocorrelation function of the stress tensor, (b) Autocorrelation function of the energy. Note the difference in 
scale of the a;-axes of the plots. 

The autocorrelation behavior of the stress tensor, however, is dramatically different from that of the energy. Note 
the difference in scale of the i-axes in Fig. ||. The autocorrelation function of the stress tensor drops so sharply 
that the exponential behavior can hardly, if at all, be seen. There is almost no sign of critical slowing down; the 
autocorrelation functions even seem to converge for larger and larger systems. Even for systems as large as 180 x 180 
spins the autocorrelation function behaves not significantly different from smaller system sizes. 

These findings can be explained as follows. The dynamic scaling hypothesis in its general form considers the 
combined spatial and timecorrelation function G(r,t), defined as 

G(r, t) = (O(r , t )O(r +rA + t)> - (O(r , t )) 2 , (63) 

for a certain operator 0(r, t). Here the dynamics of the system is explicitly taken into account by the time dependence 
of the operator. The dynamic scaling hypothesis states that 

■-2x/~t/-L-l. 



G(r,t) 



b-' x G(b~ L r, b~ z 



t), 



(64) 



where x is the critical dimension of the operator 0(r, t). In terms of this correlation function, the autocorrelation 
function g(t) of Eq. ( |59| ) can be expressed as 



9{t) 



f L d 2 TG(r,t) 



(65) 



fo <Pr G(r,0) 

The integral is over the finite volume L 2 . The dynamic scaling hypothesis (^i|) will be valid provided that the 
appearing lengths arc smaller than the correlation length £, and the times are smaller than the autocorrelation time 
r, given by r ~ £ z . For finite systems, £ ~ L. In that case, Eq. (|||) can be rephrased to 

G(r, t) = £- 2a; / z G(t- 1/z r, 1), (66) 

which yields the L- and ^-dependence in the scaling limit of Eq. (|6^). Using (|6^), the nominator is 

t (2 - 2x ^ z I d 2 rG(r,l), (67) 
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and G(r, 1) must follow the usual spatial behavior |r| 2x . Now the scaling behavior of the integral depends on whether 
it converges or diverges for large L. Making this distinction, the scaling behavior of Eq. (|65| ) is 

g(t) ~ constant for x < 1, 

g(t) ~ fori>l. (68) 

This explains our MC results: both cases indicate that g(t) must become independent of L in the scaling limit, i.e., 
for large L. The case of the energy, having x = 1, states that </(£) must converge to a value independent of t, whereas 
the case of the stress tensor implies that g(t) becomes a true power law in t. This behavior indeed can be seen in 
Fig. |J, where for the stress tensor and for the energy, the values of g(t) are plotted as a function of system size L for 
several values of t. The plot for the energy indicates that g(t) converges to 1, whereas the asymptote of g(t) for the 
stress tensor is seen to depend on t. 




FIG. 4. Plots of the autocorrelation function g(t) of Eq. ( |59| ) versus the system size L. From low to high, the plots amount 
to t = 3, t = 5 and t — 10. (a) g(t) for the stress tensor, (b) g(t) for the energy. The plots show that for the energy, g(t) 
converges to a value independent of t, which must be 1. For the stress tensor, however, g(t) converges to a value that does 
depend on t. 



The above analysis also enables us to determine the scaling of the typical simulation time of a MC simulation with 
the system size L. This scaling will depend on the MC algorithm (i.e., on z) and on the observable we want to know. 
Starting point is that we will want to obtain the same relative error in the average value (O(r)) for each system 
dimension L. The error Ao( r ) in the average is proportional to the second moment of the correlation function, 

A ° W = ^f d2^G(^,^) ■ (69) 

Here N is the total number of MC cycles, which is supposed to be much larger than the autocorrelation time L z . 
Using, as above, Eq. (|64| ) and the distinction between converging and diverging integrals with L, we obtain 

Ao(r) ~ lt LZ ~ 2X fOT X < 1 + 

\ I (70) 

A o(r) ~ JyL for x > 1 + -z. 
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The relative error is obtained by dividing these values by (O(r)), which scales as I x . The typical number of MC 
cycles N is obtained by demanding it to be such that the same relative error is obtained for all L. This yields 

N ~ L z for z > 2x - 2, 

, , (71) 
N ~ L 2x - 2 for z < 2x - 2. 



This implies that for a relevant operator, like the energy in the Ising model, faster convergence is obtained by a MC 
algorithm that has a lower value of the dynamic exponent z. However, the case of the stress tensor, x = 2, represents 
a border case, because for the Metropolis algorithm z is only slightly larger than 2. This explains why it is not 
necessary to use a more sophisticated cluster algorithm for MC simulations on the stress tensor. 

Note that the actual simulation time is, in the case of a Metropolis algorithm, proportional to L 2 N, because the 
time needed for a single MC cycle is simply proportional to the number of spins. An important consequence of the 
above is that the typical simulation time for the stress tensor is roughly proportional to L 4 . This contrasts with the 
computer time needed for transfer matrix calculations, which is exponential in L. 

Hence, in principal much larger system sizes can be reached with our MC method than in the transfer matrix 
method, to calculate the central charge. This is a promising conclusion for systems of which the value of the central 
charge up to now is an open question JlJJ]. 



VII. DISCUSSION AND CONCLUSION 



In this paper, we proposed a novel Monte Carlo technique for the calculation of the central charge and some critical 
dimensions of two-dimensional critical models. The technique is based on the universal behavior of the stress tensor, 
an operator that plays an important role in the theory of conformal invariance, but of which a lattice representation 
can easily be found as well. The rough data, following from the Monte Carlo simulation, require a decent fit program to 
extract the central charge and critical dimensions. By comparing our Monte Carlo analysis for three different models 
with their exact results, we show that the method works. We explain why, on one hand, simulations on the stress 
tensor are difficult because its expectation value is equipped with larger error bars than usual. On the other hand, it 
turns out that the simulations are much easier than usual because the stress tensor shows to be remarkably insensitive 
to critical slowing down. The latter observation notably ensures that the typical simulation time of our method scales 
with the system size L roughly as L 4 , in contrast with transfer matrix calculations, which scale exponentially as n L , 
where n is the number of different spin states. 

Hence, in principal much larger system sizes can be reached with the proposed method than with transfer matrix 
calculations. For that reason, we expect the merits of our method to lie mainly in simulations on models with a large 
number of spin states n, especially when these states become continuous, as, e.g., in the AY"-Ising model. As the 
stress tensor is highly insensitive to critical slowing down, advantages can also be obtained when no cluster algorithm 
is available for the Monte Carlo simulations. 
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